Crackling noise in three-point bending of heterogeneous materials 
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We study the crackling noise emerging during single crack propagation in a specimen under three- 
point bending conditions. Computer simulations are carried out in the framework of a discrete 
element model where the specimen is discretized in terms of convex polygons and cohesive elements 
are represented by beams. Computer simulations revealed that fracture proceeds in bursts whose 
size and waiting time distributions have a power law functional form with an exponential cutoff. 
Controlling the degree of brittleness of the sample by the amount of disorder, we obtain a scaling 
form for the characteristic quantities of crackling noise of quasi-brittle materials. Analyzing the 
spatial structure of damage we show that ahead of the crack tip a process zone is formed as a 
random sequence of broken and intact mesoscopic elements. We characterize the statistics of the 
shrinking and expanding steps of the process zone and determine the damage profile in the vicinity 
of the crack tip. 
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I. INTRODUCTION 

The brittle fracture of materials has two substantially 
different scenarios depending on the amount of struc- 
tural disorder: for homogeneous materials such as crys- 
talline solids at the critical stress a single crack is formed 
which propagates in an unstable manner. However, in 
materials with a high degree of heterogeneity fracture 
develops progressively, i.e. under an increasing external 
load first microcracks nucleate at local weaknesses which 
may then undergo several steps of growth and arrest. 
Finally macroscopic fracture occurs as the culmination 
of the gradual accumulation of damage [l|, [|J . The nu- 
cleation and growth of cracks is accompanied by the 
emission of elastic waves which can be recorded in the 
form of acoustic noise @ . Measuring acoustic emissions 
(AE) on loaded specimens is the primary source of in- 
formation on the microscopic dynamics of the fracture 
of heterogeneous brittle materials @, 0, [f| . During the 
last two decades a large number of acoustic emission ex- 
periments were carried out on different materials under 
quasi-statically increasing and constant external loads. 
These experiments revealed that the energy of acoustic 
bursts and the waiting times between consecutive events 
are characterized by power law distributions. The value 
of the exponents of crackling noise is found to be char- 
acteristic for the type of fracture, i.e. for ductile fracture 
the exponents are larger than for brittle breaking since 
large acoustic bursts are suppressed in ductile materials 

iii- 
it has long been recognized that in spite of the widely 
different length scales, acoustic emissions of fracturing 
solids and earthquakes share several common features 
3]. Recently, it has been pointed out that the probabil- 
ity distribution of waiting times T between consecutive 



earthquakes can be described by a simple scaling form 



P(T) ~ Rf(RT), 



(1) 
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where R denotes the mean rate of events in the time 
window considered. The generic law Eq. ([T]) proved to 
be valid for all geographical regions unless the time win- 
dow is sufficiently broad to fulfill the condition of sta- 
tionarity Q- Laboratory experiments on the fracture 
of different types of materials have revealed that the 
scaling form Eq. (JT|) is also valid for acoustic emissions, 
even the scaling function f(x) proved to have the same 
functional form, i.e. a gamma distribution is obtained 
f(x) = Ax 1 exp(—x / B) which means a power law decay 
followed by an exponential cutoff [8) . 

Very recently the statistical features of acoustic emis- 
sions have been analyzed during three-point bending 
tests of notched concrete specimens. The notch ensures 
that the formation of microcracks is dominated by strong 
spatial correlations in a narrow cross section of the speci- 
men. Varying the detection threshold of acoustic signals, 
for the waiting time distributions the scaling behavior 
Eq. ([IJ was recovered In Ref. [ltj the three-point 

bending experiment was modelled by discretizing the bar 
in terms of a bundle of fibers. It was a crucial feature 
of the model that after fiber breakings the load was re- 
distributed equally in the bundle, i.e. spatial correlation 
of microfractures were not taken into account. Computer 
simulation of the bending process revealed the same scal- 
ing structure Eq. ([1} of the numerical results, however, 
for the scaling function / an exponential form was ob- 
tained / ~ exp(— RT). These experimental and theo- 
retical results demonstrate the importance of the range 
of stress redistribution and correlations of consecutive 
acoustic events in the statistics of microfractures. 

Motivated by these theoretical and experimental find- 
ings, in the present paper we investigate the fracture of 
heterogeneous materials under three-point bending con- 
ditions by means of a discrete element model (DEM). 
Our two-dimensional DEM approach provides a realis- 
tic representation of the microstructure of the material, 
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the formation of microcracks, and the emerging compli- 
cated stress field naturally accounting for the correlation 
of microfractures. We analyze both the temporal evolu- 
tion and spatial structure of damage varying the degree of 
heterogeneity of the material. Our investigations showed 
that the fracture proceeds in bursts whose size and wait- 
ing time distributions have power law behavior with an 
exponential cutoff. Varying the amount of disorder of 
the sample we can control the degree of brittleness of the 
final failure of the material. A scaling form is determined 
in terms of which the distributions obtained at different 
amounts of disorder can be collapsed on a master curve. 
Ahead of the crack tip a process zone is formed composed 
of broken and intact mesoscopic elements. Our DEM ap- 
proach allows us to carry out a detailed analysis of the 
spatial structure of damage as well. 



II. DISCRETE ELEMENT MODEL FOR 
HETEROGENEOUS MATERIALS 

Three-point bending is a standard engineering test 
where a bar shaped specimen is clamped at the two ends 
and a point load is applied in the middle perpendicular 
to the longer axis of the bar. Under an increasing load 
the bar bends and finally breaks due to a crack which 
appears in the middle along the load direction. This test- 
ing method is mainly used in the engineering literature 
to characterize the quasi-static fracture strength of struc- 
tural materials such as concrete [ll| . On the other hand, 
three-point bending experiments provide an excellent op- 
portunity to study the propagation of a single crack in 
a disordered environment which is a challenging problem 
for the statistical physics of fracture Q. 

Recently, we have worked out a two-dimensional dy- 
namical model of deformable, breakable granular solids, 
which enables us to perform molecular dynamics simu- 
lation of fracture and fragmentation of solids in various 
experimental situations [12Hl5j |. Our model is an exten- 
sion of those models which are used to study the behavior 
of granular materials applying randomly shaped convex 
polygons to describe grains [12j . The initial set of poly- 
gons is obtained by the Voronoi tessellation of a rectangle 
from which specimens of appropriate shapes can be cut 
out. The average polygon size l p sets the characteris- 
tic length scale of the model system. The polygons are 
considered to be rigid bodies which can overlap when 
pressed against each other. We introduce a repulsive 
force between the overlapping particles proportional to 
the overlap area [l^ - [T5j . To capture the elastic behav- 
ior of solids we connect the unbreakable, undeformable 
polygons (grains) by elastic beams. The beams have two 
important roles in the model construction: they ensure 
cohesion and they are able to break which is essential to 
model fracture processes. The beams can be elongated, 
compressed, sheared and bent so that they exert forces 
and torques on the polygons to which they are attached. 
Figure Q] presents an example of the polygon structure 




FIG. 1. (Color online) (left) Neighboring polygons of the 
initial Voronoi tessellation are connected by beams. This way 
a triangular beam lattice is obtained, (right) Due to sub- 
sequent breaking of beams a crack forms along the edge of 
polygons. 



and the beam lattice attached to the polygons. 

In the simulations a bar shaped specimen is consid- 
ered with longer and shorter side lengths L and L Cl re- 
spectively. In order to make a realistic representation of 
three-point loading, the three loading plates are realized 
by additional polygonal elements, i.e. squares in Fig. [5] 
with side length S — 5l p much smaller than the longer 
side L = 200Z p of the bar S <C L. These loading plates in- 
teract with the particles of the bar via the overlap force, 
however, no beams are coupled to them. Strain con- 
trolled loading of the bar is implemented in such a way 
that the two loading plates at the bottom are fixed while 
the third one on the top is moved vertically downward in 
Fig. [5] with a constant speed vq. The moving plate over- 
laps the boundary polygons on the top of the bar which 
results in an increasing loading force. The stiffness of the 
plates is set high enough to keep the overlap below 20% 
of the average polygon area. Simulations were carried 
out varying the value of vo in a range, which allows for 
an efficient damping of the elastic waves and ensures a 
reasonable CPU time for the computations. The main 
advantage of three-point bending tests is that the highly 
stressed zone, where the crack appears, falls in the mid- 
dle of the bar which helps to make efficient monitoring 
of the fracture process. In order to simplify the numeri- 
cal measurements on crack propagation, we introduce a 
"weak" line in the middle of the bar in such a way that 
solely those beams are allowed to break which connect 
the two sides of the line (see Fig. 



A. Disordered beam breaking 

The beams, modeling cohesive forces between grains, 
can be broken according to a physical breaking rule, 
which takes into account the stretching and bending of 
the connections [H, H3] 
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Here e denotes the longitudinal deformation of a beam, 
while Oi and O2 are bending angles at the two beam 
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FIG. 2. (Color online) Three-point bending of a bar com- 
posed of polygonal particles. The particles are coupled by 
elastic beams which are colored according to the longitudinal 
deformation (yellow: nearly unstressed beams; red and black: 
elongated beams; blue and green: compressed beams). Beams 
are allowed to break solely along the center line of the bar. A 
relatively small sample is presented to have a clear view on 
the details of the model construction. The two loading plates 
at the bottom are fixed while the third one on the top moves 
downward. 



ends. The breaking rule Eq. contains two parameters 
£th, @th controlling the relative importance of the stretch- 
ing and bending breaking modes, respectively. The en- 
ergy stored in a beam just before breaking is released in 
the breakage giving rise to energy dissipation. At the bro- 
ken beams along the surface of the polygons cracks are 
generated inside the solid and as a result of the successive 
beam breaking the solid falls apart (see Fig.[T]). The time 
evolution of the polygonal solid is obtained by solving the 
equations of motion of the individual polygons. At each 
iteration step we evaluate the breaking criterion Eq. ^ 
and remove those beams which fulfill the condition. The 
simulation is continued until all beams break along the 
weak interface and the specimen falls apart. (For more 
details of the model construction see Refs. fl2l. Il3j.) 

The breaking parameters Eth and Qth of beams are 
stochastic variables in the model, i.e. they are sampled 
from probability density functions p(sth) and p(Qth)- 
The Weibull distribution provides a comprehensive de- 
scription of the stochastic fracture strength of materials, 
hence, for both threshold values the Weibull form is pre- 
scribed 



range 1 < m < 50 for both threshold distributions in or- 
der to control the amount of disorder in the system. We 
note that the elastic constants of beam elements depend 
on the Young's modulus of beams, furthermore, also on 
their length and cross section. In the model the geome- 
try of beams is determined by the Voronoi tessellation, 
i.e. the length and cross section of beams are defined as 
the distance between the center of mass and the length of 
the common side of the two neighboring polygons, respec- 
tively. It has the consequence that besides the strength 
disorder of beams there is also structural disorder in the 
system determined by the initial Voronoi tessellation. 



III. MACROSCOPIC RESPONSE 

In our strain controlled three-point bending experi- 
ment the mechanical response of the material can be 
characterized numerically by measuring the force F act- 
ing on the moving plate at the top of the sample as a 
function of time t (see Fig. [5]). Figure [3] presents the 
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where x denotes the two breaking thresholds E t h,®th- 
The Weibull distribution has two parameters: A sets the 
characteristic scale of threshold values while the expo- 
nent m determines the scatter of the variable. Increasing 
the value of the exponent m the width of the Weibull 
distribution Eq. © decreases and converges to the delta 
function in the limit m — > oo. Varying the scale param- 
eters A £ and Ae of the breaking thresholds the relative 
importance of stretching and bending can be controlled in 
the breaking process. For clarity in the present paper we 
only investigate the two limiting cases of beam breaking 
dominated by pure stretching or bending with the param- 
eter settings A e = 0.05, Ae = 100, or A 6 = 100, Ae = 1, 
respectively. The Weibull exponents were changed in the 



FIG. 3. (Color online) Force F normalized by the maximum 
value Fmax as a function of time t during the loading process. 
Since the deflection of the bar is proportional to t the curve 
can be considered as the constitutive curve of the sample. Os- 
cillations occur due to elastic waves generated by the loading 
process. After the peak the decreasing part of F(t) indicates 
stable crack propagation where our analysis is focused, ttot 
denotes the time of the last beam breaking. 



force-time curve obtained for a sample with the Weibull 
exponent m = 2.0 in the stretching limit. Since the load- 
ing plate moves at a constant speed the deflection of the 
bar is proportional to t so that F(t) can be considered to 
be the constitutive curve of the sample. It can be seen 
that the macroscopic response is linear all the way up 
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to the peak, where the force drops suddenly. This drop 
becomes more and more drastic as we increase the brittle- 
ness of the sample by increasing the value of the Weibull 
exponent. At the beginning of the loading process, the 
smooth oscillations about the linear in the constitutive 
curve arise due to elastic waves generated by the loading 
plate. As the force approaches its maximum, the curve 
becomes more and more noisy due to internal damage 
being accumulated in the form of microcracks nucleating 
throughout the breakable interface. It can be observed 
in Fig. [3] that after the maximum, the force drops rather 
drastically, however, the failure is not totally abrupt. The 
sharp drop-down of the force is followed by stable crack 
propagation where the crack gradually advances until the 
sample falls apart. 

As the bar is loaded at a constant strain-rate, micro- 
cracks - corresponding to uncorrelated beam breaking 
events - start nucleating throughout the interface. This 
way damage is accumulated inside the sample before the 
onset time of crack propagation. Local beam breaking 
inside the sample is always brittle, however, the disor- 
der of breaking thresholds can result in a quasi-brittle 
macroscopic response where the constitutive curve ex- 
hibits a non-linear behavior. Due to the disturbing ef- 
fect of elastic waves, numerically it is difficult to quantify 
the strength of non-linearity of the F(t) curve in Fig. [3] 
Hence, we characterize the degree of brittleness of the 
sample and its dependence on the amount of threshold 
disorder by measuring the accumulated damage prior to 
the peak of the force as a function of the Weibull expo- 
nent m. Figure @] shows the damage parameter D of the 
model, defined as the fraction of beams broken before 
a single crack starts propagating, for the stretching and 
bending limits, as a function of the Weibull exponent m. 
It can be observed that the curves can be very well fitted 
with the functional form 

D(m) = B + Am- fi , (4) 

where all parameters A, B and /_i proved to be different 
for the stretching and bending limits. The power law 
form of D(m) can be motivated by the following simpli- 
fied assumption: Let us consider a mean-field approxima- 
tion of the system, where all the beams along the inter- 
face share the same e strain at any given time during the 
process. This way the breakable interface of the bar is 
substituted by a parallel bundle of beams with equal load 
sharing, whose breaking process can easily be described 
analytically [l6l - fl8j . The fraction of intact beams at any 
e can be given as 1 — P{e) where P(s) is the cumulative 
probability distribution of the breaking thresholds. The 
macroscopic stress a as a function of strain e can then 
be written in the form 

tr(e) - [1 - P(e)] Ee = e -( £ / A ) m Ee, (5) 

where E is the Young's modulus of the beams. Under 
strain controlled loading of the bundle, stable crack prop- 
agation starts at the peak of the constitutive curve cr(e). 
After differentiating Eq. ([5]) the position of the maximum 




FIG. 4. Damage accumulated up to the peak of F(t) as a 
function of the Weibull exponent m. The curves can be very 
well fitted with the functional form Eq. @. The value of 
the exponent is fi = 1.9 and /i = 1.5 for the stretching and 
bending limits, respectively. 



e c reads as e c = A(l/m) 1 /™ 1 for the Weibull distribution. 
The fraction of broken beams accumulated up to the peak 
of cr(e) can be obtained by plugging e c into the cumula- 
tive distribution of thresholds P(s c ), hence, the damage 
parameter D as a function of the Weibull exponent m 
can be cast into the final form for large enough m values 

D(m) w l-e^ 1 /™) ~m _1 . (6) 

The numerical results on the amount of damage prior to 
the force peak in Fig. 2] are consistent with the above 
analytic prediction. The higher value of the measured 
exponents ^ s = 1.9 ± 0.1 (stretching) and fj, b = 1.5 ± 
0.1 (bending) is the consequence of the strain gradient 
in the load direction, which was completely neglected in 
the analytic calculations. Note that in the limit of high 
to values, the amount of damage does not converge to 
zero, instead it takes a finite value B > 0. The non-zero 
value of B in Eq. ^ can be attributed to the structural 
disorder in the sample, which is present and is the same 
for all values of the Weibull exponent. This structural 
disorder gives rise to fluctuations of the beam parameters 
which in turn result in a noisy breaking sequence in spite 
of the constant breaking parameters [1J, [l4| . 

Perfectly brittle failure of the bar would be charac- 
terized by a linear behavior of F(t) up to the maxi- 
mum without any damaging which is then followed by 
an abrupt breaking. Our simulation results demonstrate 
that varying the amount of threshold disorder we can 
control the degree of brittleness of the DEM sample from 
highly (but not perfectly) brittle to quasi-brittle. It is a 
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FIG. 5. Time series of bursts in a single fracture simulation. 
The bursts are correlated breaking sequences of beams which 
then result in sudden jumps of the extending crack. N denotes 
the total number of beams along the weak interface where the 
crack propagates. For all the simulations its value was set to 
N — 200. At the beginning of the loading process, for a 
considerable time no breaking occurs, most of the breaking 
events appear at larger deflections beyond the peak of the 
constitutive curve (see Fig. [3}. Hence, we magnify the final 
section of the bending process. 



very interesting question how the degree of brittleness 
affects the properties of crackling noise and the spatial 
structure of damage along the interface. 



IV. CRACKLING NOISE DURING CRACK 
PROPAGATION 

A snapshot of the computer simulation of a three-point 
bending test is presented in Figure [SJ The color code 
shows that the bottom of the specimen is highly elon- 
gated that's why the crack starts here. In the vicinity of 
the crack tip the beams are strongly elongated indicat- 
ing a high stress concentration ahead of the crack which 
provides the driving force for crack propagation. At the 
top of the bar the color code indicates that compressive 
stresses arise. The constant speed of the loading plate 
ensures a strain controlled loading of the specimen at a 
fixed strain rate. The low value of the loading speed has 
the consequence that in each iteration step of the molec- 
ular dynamics simulation either no beam breaking occurs 
or only a single beam breaks. After a local breaking event 
the stress gets redistributed increasing the stress concen- 
tration on the intact elements ahead of the crack. This 
load redistribution may give rise to additional breakings 
resulting in a correlated trail of breaking events. 

In order to identify bursts of local breakings we in- 



FIG. 6. Inset: Avalanche size distributions for the absolute 
stretching limit varying the value of the Weibull exponent m. 
The main panel presents the excellent data collapse obtained 
by rescaling the distributions with the average burst size ac- 
cording to Eq. Scaling exponents: a s A = 1.4 ± 0.5, /J^ = 
1.8 ± 1. The parameter values obtained by fitting Eq. ((8]) are 
a s A = 0.55, t! = 1.3 ± 0.2, b s A = 2.2, 8%. = 1.5 ± 0.3. 



troduce a correlation time t corr : if the time difference of 
two consecutive beam breakings occurring at times U and 
tj_l_i is smaller than the correlation time t i+ i — tj < t corr 
the two breakings are considered to belong to the same 
burst. The value of the correlation time was chosen in 
such a way that it is larger than the time step At used in 
the integration of the equation of motion but it is much 
smaller than the total duration t tot of the breaking pro- 
cess, i.e. we set t corr = 10 At for which 10°t corr < t to t 
holds. The size of bursts A is defined as the num- 
ber of beams breaking during the correlated sequence. 
When the amount of disorder is very high m — > 1, es- 
pecially in the bending limit of breakings, it may hap- 
pen in DEM simulations that very distant beams break 
within the correlation time, however, without any corre- 
lation. To obtain information on the strength of spatial 
correlations in an avalanche, we calculate the distance 
hj = \Uj ~~ between consecutive beam breakings 

with the positions yj and yj+i and sum it up inside 
an avalanche ft — ^2f = i hj. For a strongly correlated 
avalanche where each consecutive breaking occurs on ad- 
jacent beams the ratio of h and of the burst size A is 
close to the characteristic polygone size ft/ A « l p . In 
order to filter out avalanches dominated by random coin- 
cidences we introduce a threshold value for this ratio, i.e. 
those avalanches for which ft/A > 21 p holds are removed 
from the statistics. Computer simulations showed that 
in the stretching limit the above condition has no effect, 
however, in the bending limit where a high amount of 
distributed cracking occurs, about 10% of the avalanches 
are filtered out due to random coincidences (compare also 
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FIG. 7. Inset: avalanche size distributions for the absolute 
bending limit. The main panel shows that rescaling the dis- 
tributions according to Eq. J?} an excellent data collapse is 
obtained. Scaling exponents: a b A — 1.4 ± 0.5, /3 A = 1.8 ± 1. 
The fit parameters of the scaling function are a b A = 0.85, t a = 
0.8 ± 0.3, b\ = 1.4, 5 b A = 1.3 ± 0.3. 



to Fig. gl). 

Figure [5] presents the size of bursts in a single frac- 
ture simulation at the time of their appearance. It can 
be seen in the figure that the bursts are separated by 
silent periods with variable length. These waiting times 
between bursts characterize the duration of states where 
the crack tip is pinned due to the presence of some strong 
beams. At the beginning of the loading process the bursts 
are small compared to the cross section of the specimen 
(maximum crack length), however, with increasing de- 
flection of the bar the burst size A increases and reaches 
a maximum somewhat before the last breaking. After 
the maximum the burst size decreases showing that as 
the crack approaches the top of the bar it slows down 
due to the compressive zone. 

We determined numerically the size distribution of 
bursts P( A) varying the amount of disorder in the failure 
thresholds. The size distribution obtained at different 
values of the Weibull exponent is presented in the in- 
sets of Fig. [5] and Fig. [7] for the absolute stretching and 
bending limits, respectively. It is interesting to note that 
increasing the Weibull exponent to, i.e. decreasing the 
amount of disorder, the bursts get larger but the func- 
tional form does not change. For small bursts a power 
law behavior is obtained followed by a rapidly decreas- 
ing cutoff regime. The main panels of Fig. [6] and Fig. [7] 
demonstrate that using the average burst size A as a scal- 
ing variable, the burst size distributions -P(A) obtained 
at different m values can be collapsed on a master curve. 



The data collapse implies the scaling structure 

P(A) = A^(A/A Q ), (7) 

where the values of the exponents were determined nu- 
merically a% = 1.4±0.5, A = 1.8±1, and a\ = 1.4±0.5, 
j3 b A = 1.8 ± 1 which provide the best quality collapse for 
stretching and bending, respectively. The scaling func- 
tion / can be very well fitted by the form 

f(x)=ax- T e-( x / b *>\ (8) 

where the parameter values providing the best fit are 
a s A = 0.55, t a = 1.3 ±0.2,6^ = 2.2, S S A = 1.5 ± 0.3 
(stretching), and a b A = 0.85, t a = 0.8 ±0.1,6^ = 
1A,S A = 1.3 ± 0.3 (bending). The results demonstrate 
that the growth of the crack is not a smooth process, the 
slow driving results in a jerky crack propagation which is 
composed of a large number of discrete steps. The growth 
steps are sudden outbreaks with a variable length. The 
correlation of consecutive local breakings leads to a power 
law functional form limited by an exponential cutoff. The 
most interesting outcome of the calculations is that the 
amount of disorder only affects the characteristic scale 
of bursts but the functional form and the value of the 
power law exponents t a and t a remains the same. We 
note the differences between the limits of stretching and 
bending dominated breaking, especially the deviation of 
the exponents t a and t a is beyond the error bars. The 
functional form Eq. (|8|) has also been found to provide 
a good quality description of the amplitude distribution 
of acoustic bursts in three-point bending experiments on 
concrete samples 0, H3] • 

It can be observed in Fig. [5] that the bursts are sep- 
arated by silent periods where no beam breaking oc- 
curs. The advancing loading plate gradually increases 
the load on the system which reactivates the crack after 
some waiting time T. It can be seen in Fig. [5] that the 
duration T of these waiting times can vary in a broad 
range. In Figure [5] the waiting time distributions P{T) 
are presented for the stretching limit separated for high 
(inset) and low disorder (main panel). It is interesting 
to note that for low enough disorder (main panel of Fig. 
[8]) the distributions P(T) are all the same, no depen- 
dence on the Weibull exponent could be pointed out. 
The functional form of P(T) can be very well fitted by 
the expression Eq. ([5]) where the value of the exponent 
Tj. = 1.9 ± 0.15 was obtained. The relatively high value 
of implies that long waiting times are very rare in the 
trail of bursts when the material is very brittle. However, 
in the limit of high disorder m — » 1 (see inset of Fig. [5]) 
waiting times span a broader range and reach an order 
of magnitude larger values than for the very brittle ma- 
terials with low disorder. The most remarkable feature 
of waiting time distributions is that increasing the disor- 
der the exponent of the power law regime changes to the 
lower value t|> = 1.5 coinciding with the recurrence time 
exponent of one-dimensional random walks. In the ab- 
solute bending limit (see Fig. [SJ P{T) has qualitatively 
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FIG. 10. (Color online) Crack tip and the process zone in 
front of it. Red dots indicate the positions of broken beams. 
The process zone is identified as a sequence of broken and 
intact beams starting at the crack tip and ending at the start 
of the continuous sequence of intact fibers. 



V. SPATIAL STRUCTURE OF DAMAGE 



FIG. 8. Waiting time distributions for absolute stretching. 
The main panel presents curves for low disorder, where the 
fit was obtained with the exponent r-f = 1.9 ±0.15. The 
inset shows the corresponding curves for high disorder, where 
a crossover is obtained to a lower exponent r-f = 1.5 ± 0.1. 




T/ttot 

FIG. 9. Waiting time distributions for the absolute bending 
limit. The amount of disorder only affects the cutoff but the 
exponent is constant r|> = 1.8 ± 0.15. 



It has been discussed in Sec. lIIII that at the time of the 
force drop in Fig. [3J a crack initiates along the breakable 
interface of the sample and proceeds in a jerky way. The 
crackling noise analyzed in the previous section charac- 
terizes the temporal fluctuations of the advancing crack. 
A very important advantage of our modeling approach 
is that it allows us to investigate the spatial structure 
of damage as well. Based on the sample geometry and 
loading conditions illustrated in Fig. f2J the crack can be 
identified as a continuous region of broken beams start- 
ing from the bottom of the interface. The high stress 
concentration ahead of the crack tip and the quenched 
disorder of the local strength of beams give rise to the 
emergence of a sequence of broken and intact beams fol- 
lowed by a continuous region of intact elements. Figure 
[TUl illustrates that in the framework of our discrete ele- 
ment model the crack tip can be precisely defined as the 
position of the first intact beam starting from the bot- 
tom of the specimen. The sparse region of broken and 
intact beams, extending from the crack tip to the last 
broken fiber, can be identified as the fracture process 
zone (FPZ) whose dynamics has a strong influence on 
the time evolution of the breaking process . We note 
that the extension of the process zone is also affected by 
the background damage D nucleated in an uncorrelated 
manner before the crack starts (see Fig. 0]). The higher 
amount of background damage in the bending limit re- 
sults in a larger extension of the process zone than for 
the case of stretching dominated breaking. 



the same behavior as in the stretching limit. Due to the 
fragility of the system at all Weibull exponents to, the 
change of disorder only results in a change of the cut- 
off, however, the value of the exponent of the power law 
regime is constant Tj, = 1.8 ± 0.15. 



A. Dynamics of the process zone 

The dynamics of the process zone strongly determines 
the advancement of the crack. As the beams break, the 
extension of the process zone changes in discrete steps: 
when a new breaking nucleates inside the intact zone the 
FPZ extends by a length l nuc i called nucleation length. 
The process zone shrinks when the beam at the crack tip 
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breaks resulting in a jump of the crack tip (CTJ) with 
a distance Ictj- I n order to characterize the dynam- 




10" 2 10" 1 

Ictj/ L c 

FIG. 11. Crack tip jump length distributions for the absolute 
bending limit. The values of the exponent j b of the fitted 
curves are 1.8 and 4. 

ics of FPZ we investigate the probability distribution of 
the nucleation length p{l nuc i) and of the length of crack 
tip jumps pQctj)- We note that the notion of crack 
tip jump length has also recently been introduced in the 
framework of Quantized Fracture Mechanics j20l - l22| . It 
can be observed in Fig. [TT] that for low disorder (high 
Weibull exponent m) in the bending limit of breakings 
the distribution of crack tip jumps has a power law decay 
in the regime of small Ictj values which is complemented 
by an exponential form for large Ictj 

p(l CTJ ) = al^ TJ + ce~ lcT >/ b . (9) 

The additive coupling of the two terms of Eq. shows 
that the small and large crack tip jumps are generated 
by different mechanisms, i.e. the small ones are deter- 
mined by the stress concentration at the crack tip and 
by the resulting correlation of local breakings. However, 
the Poisson-like behavior of very large crack tip jumps 
captured by the second term of the equation originates 
from the randomness of the initial jump-in of the crack 
at the onset of crack propagation. The value of the ex- 
ponent 7 & changes from j b = 1.8 ± 0.1 to j b = 4 ± 0.2 as 
the amount of disorder decreases. In the stretching limit 
of breakings the sample behaves in a less fragile way ac- 
cumulating less background damage before the onset of 
crack propagation (see also Fig. Hence, in Fig. Q2] 
we obtain a power law decay with an exponential cutoff 
but the additive exponential term does not occur. Sim- 
ulations showed that the value of the exponent does not 
depend on the amount of disorder 7 s = 2.2 ± 0.1 but 
it falls between the low and high disorder limits of the 
bending counterpart. 




Ictj/ L c 

FIG. 12. Crack tip jump length distributions for the absolute 
stretching limit. The value of the exponent of the power law 
regime is 7 s = 2.2 ±0.1, it does not depend on the amount of 
disorder. 



The statistics of crack tip jumps was accumulated over 
the entire time evolution of the bending process. In order 
to obtain information about the nucleation length l nuc i, 
we analyze microcracks occurring ahead of the crack tip 
at the time when the crack spans approximately half of 
the cross section of the specimen. It can be observed in 
Fig.[l3]that in the bending limit the distribution p(l nuc i) 
hardly changes when the amount of disorder is varied. A 
power law form is obtained for low length values 

P(lnucl) ~ Cd (10) 

with an exponent K b = 1.8 ± 0.15. When beam breaking 
is dominated by tensile stresses the situation drastically 
changes: at high disorder the nucleation distance can ex- 
tend up to 40% of the cross section of the specimen as 
in the bending limit, i.e. very remote beams may also 
break if they are weak enough. At low disorder, how- 
ever, the nucleation of new broken beams gets localized 
to the close vicinity of the crack tip where the maximum 
of lnuci extends only up to 1 — 2% of the cross section. 
These results on the distance to new nucleations and on 
the length of crack tip jumps clearly show that in the 
case of bending dominated breaking varying the amount 
of disorder does not have a strong effect on the spatial dis- 
tribution of damage. The dynamics of the process zone 
is mainly determined by the long range redistribution of 
stresses arising from the bending distorsion. However, 
when breaking is dominated by tensile deformation, dis- 
order plays a crucial role in the evolution of the fracture 
process zone, i.e. at low disorder the process zone expands 
in a large number of small steps while shrinking occurs 
in the form of a few larger jumps. When the disorder 
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is high, both shrinking and expanding steps can span a 
broad range. 




nucl. 



For the exponent oj linear fracture mechanics predicts 
the value oj = 1/2 [l9j], while fractal cracks and plastic 
or hyper-elastic constitutive laws lead to different values 
of oj [23[ . Hence, the probability of beam breaking as a 
function of r can be estimated as 

d(r) - P(a(r)), (13) 

where P(x) is the cumulative distribution function of the 
breaking thresholds. Since our breaking thresholds are 
Weibull distributed, wehaved(r) = l- e -( <7 ( r )/ A ) m = 1- 
e _ ( br ) where b = (a/ A)" 1 . Restricting the calculation 
for small distance we arrive at the form 

d(r) « br- um . (14) 

The curves of Eq. (fTT|) in Fig. [15] are consistent with 
the analytical expression Eq. (fT4"|) for low values of r. 
Comparing the results to Eq. (fTTj) the exponent p ob- 



FIG. 13. Nucleation length distributions for the absolute 
bending limit varying the Weibull exponent m. The amount 
of disorder does not have a relevant effect on the distribution. 
The exponent of the fitted power law is k 6 = 1.8. 



B. Damage profile 

To obtain quantitative estimates for the extension of 
the process zone, we calculated the spatial distribution 
of damage in front of the crack tip when the crack spans 
half of the specimen's cross section. Fig. [15] presents 
the damage d, i.e. the probability of beam breaking as 
a function of the distance r measured from the crack tip 
for the stretching limit of the model. The numerical re- 
sults clearly demonstrate that larger Weibull exponents, 
i.e. higher degree of brittleness results in smaller process 
zones. For all Weibull exponents the curves can be well 
described by a power law functional form with an expo- 
nential cutoff 



d(r) ~ r p exp (r/r ), 



(11) 



where the extension of the process zone can be charac- 
terized by the length tq. It can be observed in Fig. [15] 
that both the exponent p and the characteristic length ro 
depend on the amount of disorder m. Fitting the formula 
Eq. (fTTj) to the simulated data we obtained the following 
parameter values: p — 0.5, r — 0.05 (m = 2), p = 1.0, 
r = 0.035 (m = 4), p = 1.5, r = 0.033 (m = 6). 

In order to obtain an analytic understanding of the 
damage profile d(r) we can start from the result of frac- 
ture mechanics that in the vicinity of the crack tip the 
stress has a power law decay 




cr(r) 



(12) 



FIG. 14. Nucleation length distributions for the absolute 
stretching limit. Decreasing disorder leads to localization of 
damage at the crack tip. 



tained by fitting the numerical data can be written as 
a product of the exponents of stress decay and disorder 
p = ojm. Substituting the numerical values of p and the 
Weibull exponents m the exponent oj describing the de- 
cay of the stress field can be determined as oj « 0.25 for 
all to values. It is important to emphasize that this value 
of oj falls quite close to the analytic result of oj = 1/2 of 
linear fracture mechanics [l9j], furthermore, the indepen- 
dence of the numerically obtained oj from the Weibull 
exponent to shows the consistency of the results. The 
remarkable feature of the results is that the amount of 
disorder can have a strong effect both on the shape and 
extension of the process zone. 
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VI. DISCUSSION 

The emergence of crackling noise is a ubiquitous phe- 
nomenon of the fracture of heterogeneous materials which 
can also be exploited to monitor the time evolution of the 
fracture process. Theoretical studies of crackling noise 
are usually based on stochastic fracture models such as 
fiber bundles HHHH and the fuse model 0. As a novel 
approach to the problem, in the present paper we inves- 
tigated the properties of crackling noise emerging during 
the jerky propagation of a crack in three-point bending 
tests using a discrete element modeling technique. Two 
limiting cases of breakings were analyzed where stretch- 
ing or bending dominates the local breaking of beams. 
We proposed a numerical technique to identify avalanches 
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FIG. 15. Damage profile for different values of the Weibull 
exponent m. The distance r measured from the crack tip 
is normalized by the cross section L c of the specimen. The 
values of the fitting parameters are: p — 0.5, ro = 0.05 (m = 
2), p = 1.0, r = 0.035 (m = 4), p = 1.5, r = 0.033 (m = 6). 



concrete specimens d, [To|, [H|-[3(| . This agreement also 
demonstrates the importance of spatial correlations of 
consecutive microfractures in the emergence of crackling 
noise. 

An important advantage of our DEM approach is that 
it provides direct access to the spatial structure of dam- 
age. Simulations revealed that ahead of the crack tip a 
process zone develops which is a sparse region of broken 
and intact elements. The fracture process zone proved to 
play an important role in the advancement of the crack: 
on the one hand the crack progresses by shrinking and 
expanding steps of the zone, on the other hand, micro- 
cracks can shield the stress field around the crack tip 
which helps to stabilize the system. Recently the spatial 
structure of damage has been analyzed in the framework 
of the fuse model [3l|, [32|. Quasistatic loading simula- 
tions were performed starting with a notch in the middle 
of the fuse lattice analyzing the damage structure in the 
vicinity of the crack tip just before macroscopic break- 
down. It was found that the damage profile has an expo- 
nential decay along the line of the crack and the charac- 
teristic length scale of the exponential was suggested as 
the extension of the process zone. Since linear fracture 
mechanics predicts a power law decay of the stress to 
the background level ahead of the crack tip, the authors 
argued that the cloud of microcracks shields the crack 
tip giving rise to an exponential decay. In our system at 
short distances a power law decay of the damage profile 
was obtained which is followed by an exponential cutoff. 
We think the power law functional form prevails in our 
system for the damage profile because microcrack nucle- 
ation cannot occur in the two-dimensional plane but it 
is restricted to a weak "line" in the sample which de- 
creases the effect of shielding. This shielding, however, 
is responsible for the lower exponent of the stress decay 
lo and for the exponential cutoff of the damage profile. 
Computer simulation are complemented by analytic cal- 
culations under simplifying conditions, which provided a 
reasonable agreement with the numerical results. 



based on the temporal and spatial correlation of micro- 
fractures. 

We showed that for quasi-brittle materials the size of 
bursts and the waiting times between consecutive events 
are characterized by power law functional forms with an 
exponential cutoff. The numerical value of the exponents 
have a reasonable agreement with recent experimental 
findings on crackling noise in three-point bending tests on 
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